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Abstract 



Boundary conditions strongly affect the results of numerical compu- 
tations for finite size inhomogeneous or incommensurate structures. 
We present a method which allows to deal with this problem, both 
for ground state and for critical properties: it combines fluctuating 
boundary conditions and specific histogram techniques. Our approach 
concerns classical systems possessing a continuous symmetry as well as 
quantum systems. In particular, current-current correlation functions, 
which probe large scale coherence of the states, can be accurately eval- 
uated. We illustrate our method on a frustrated two dimensional XY 
model. 
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Owing to their speed and power, computers are nowadays considered as experi- 
mental tools in condensed matter physics. One of their limitations stems from the 
fact that they can only model finite size systems and it is a drawback in some cases: 
for instance, in the study of critical phenomena one needs to take the infinite size 
limit. To overcome this problem one may specify boundary conditions (BC). One 
of the most common BC consist in imposing periodic boundary conditions (PBC). 
Unfortunately PBC may lead to numerical errors and to frustration effects. For 
classical hamiltonians on a lattice this is the case when incommensurate ground 
states are favored or when interactions are long rangeBi. Dipolar magnets fall into 
that latter class: surface contributions to the energy of the system lead to inho- 
mogeneous structured!!. For quantum systems, current carrying states are affected 
by boundaries (e.g in mesoscopic ringsi). The purpose of this work is to present 
a method designed to handle numerically systems possessing a continuous sym- 
metry and sensitive to boundary effects. It combines self-determined (fluctuating) 
boundary conditions (FBC)iS and specific histogram techniques. The latter feature 
allows to study both the ground state and the critical regime of inhomogeneous or 
incommensurate structures. In particular current-current correlation functions are 
obtained straightforwardly and the method lends itself easily to correction to scaling 
analysis in the critical regime. For the sake of simplicity, the main characteristics of 
the boundary condition histogram technique are presented for classical XY spins on 
a D dimensional lattice of linear size L. The validity of our approach goes beyond 
XY systems, and it concerns the quantum case as well. 

We begin with an analytic derivation of the method; starting with the fluctu- 
ating boundary ensemble, we introduce observables. This approach is the basis of 
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a numerical study of incommensurate phases, including the case when the modula- 
tion varies with external parameters. However, near commensurate-incommensurate 
transitions or in the case of inhomogeneous ground states, large numerical errors 
occur. We then introduce a specific histogram technique (appropriate for the fluc- 
tuating boundary ensemble) and we show how it can be used to probe the nature of 
equilibrium states and to gain access to thermodynamic quantities such as current- 
current correlation functions, including in the critical regime. Next we illustrate the 
features of the boundary histogram technique on the row mo 

dellH 

a frustrated 

anisotropic 2D XY model: indeed, for some anisotropy range, the state of the sys- 
tem evolves from an incommensurate state to an inhomogeneous phase upon raising 
the temperature. 



I. ANALYTIC APPROACH 

We consider two components {XY) spins S(i) on a D dimensional lattice of 
linear size L. PBC impose an L-periodicity to the system: 

S(ii, + L, ...,i D ) = S(ii, ...,i D ) ; V(ii, ...,i D ) (1) 

As mentioned earlier, when the spins want to form an incommensurate structure, 
PBC generate frustration and introduce systematic numerical errors. Also, inasmuch 
as the pitch of the spiral may vary with T and other parameters, these BC have to 
smoothly evolve when one changes external parameters like the anisotropics or the 
temperature^. 

(Self-consistent) Fluctuating Boundary Conditions have been proposed to over- 
come these problems!'!. The main feature of FBC is to add new dynamical variables 
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A a (a = 1,2, ...,D ) corresponding to a shift at the boundaries. In equilibrium the 
new "boundary variables" A a will fluctuate around their most probable value A° . 

Variants of the FBC method have also been used to accelerate the approach to 
the asymptotic regime, by removing some unwanted correction to scalingS. This 
improvement applies both to ferromagnetic and non-ferromagnetic systems. 



A. Fluctuating boundary conditions 

The partition function of an L D system of XY spins with PBC is: 

Z PBC = j ... llX[i*,eA-^<^--^) (2) 

where the PBC are expressed by the constraint 

$(f + n x Lu x + n 2 Lu 2 + ... + n D Lu D ) = $(f) (3) 

for any set of integers ni,n 2 , ...,n D . This condition can be fulfilled by putting the 
lattice on a torus. 

The FBC method imposes instead the following constraint 
$(r + n x Lux + n 2 Lu 2 + ... + n D Lu D ) = $(f) + riiLAx + n 2 LA 2 + ... + n D LA D 

(4) 

where A 1; A^> are new dynamical degrees of freedom, corresponding to a shift at 
the boundaries. Note that using FBC allows us to preserve translational invariance 
(contrary to the free BC). Performing a change of variables 

$(f) = cp(r) + A • f (5) 
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with A = Ai^i + ... + AdUd, the constraint on ip becomes 

(p(r + n 1 Lu 1 + ...+n D Lu D )=(p(f) (6) 

In terms of the new variable ip the partition function of the L D system with FBC 
is: 




Wrapping the lattice around a torus automatically enforces the constraint (Eq||). 
It is important to note that the integration in Eq.|7| is over an interval of size 27r for 
(fi , whereas it is over a 2ir/L interval for Ai, Ad- The integration range reflects 
the periodicity of the Hamiltonian. 

Z FBC can be factorized as a product of a set of partition functions, Z(A), each 
one corresponding to a fixed shift A at the boundaries. 

Z FBC = L D ^ Z(A)d D A = L D ^ d D Ae-? LD ^ (8) 

J-it/L J-n/L 

where /(A) is the ^ periodic free energy density associated with the shift A at the 
boundary: 

f(A) = -T\n(Z(A))/L D (9) 

For A = 0, we recover the PBC case , i.e. Zp^c — Z(A = 0) and 
fp BC = /(A = 0). 

1. The low temperature regime 

First let us consider the case of a "ferromagnetic phase" at low temperature. 
The integral (Eq||) is dominated by the minima of /(A); so here, A = . We are 
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then allowed to make an expansion of the integrand around A = 0; for a lattice with 
inversion symmetry, we get: 



J FBC 



L D 



-tt/L 



d D A e 



A?- 



1 <S 2 /(A) 

2 SAl 



D 



(10) 



The second derivatives of the free energy are related to the components 7 11 , ...,7 



DD 



of the spin rigidity tensor 7 by a geometrical factor p 



5A( 



5A 2 



D 



(e.g, in 2D, p is 1 for the square lattice, and ^= for the triangular lattice). 
We rewrite Eq.llO] , after rescaling A, A' = LA : 



Z FBC ~Z PBC I d D A'exp 



-Y/L D ~ 2 (l 11 A' 2 1 + ... + 1 DD A'l 



(12) 



We deduce the formal expression connecting Zfbc an d Zpbc at low temperature, 
by extending the domain of integration of A' to [—00, +00] : 



jfbc — Apse 



^ + Q(J_ 



(13) 



,,DD 



where 7 = ^7 Xi • 7" • ... ■ 7 X 

For a system with a helical phase at low temperature, the same analysis can be 
repeated, leading to the expressions: 



J FBC 



Z(A C 



d D A'exp 



(14) 



and 



A 



5A 2 D 



(16) 



A° 



Here A is the minimum of /(A) , determined modulo -j- from: 



<5Ax 



*/(A) 



A° 



<5A 



(17) 



A° 



For a spiral phase, the pitch Qo is the 4r determination of A such that <£>(r) ~ 
in equilibrium (see Eq.[|). 
From Eq. [14] we deduce! 
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li 



where Va = 1, ..,D xa„ = P < (A Q — A°) 2 > is the susceptibility for A a . The 
validity of Eq.|18| is restricted to the low temperature regime; more precisely, to the 



DD 



P 



18) 



region of the phase diagram where Vet 



1,..,D p~f aa »l. 



2. Observables: 

Using FBC implies that thermodynamic quantities such as the energy or the 
spin rigidity will be different from those evaluated with PBC From our previous 
discussion, this difference shows to leading order in size in the incommensurate 
phase but it is also present for a commensurate phase (to order 0(1/L D )). 

For the sake of simplicity we assume in this section that /(A) has its minimum 
at A = . The results, however, do not depend on the particular value of A which 
minimizes /. 

For a given observable 0($j) the values computed with PBC and with FBC are: 
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(0) PBC = -^—1.1 Hd^O(^)e^ Jijm< ^ j} (19) 

^PBC J J —n j 

( ) FBC = — ^ — L D [ n/L d D A [ ... f n^0(^- + A •f.Oe^^ J - cosfe -^ +A - (rl ^ )) 

ZfBC J-n/L J J-n i 

(20) 

The two expressions generally differ by an 0(1/L D ) term: for instance the aver- 
age energy density is obtained by taking the derivative of Eq.|l5] with respect to (3 
yielding 



e F BC = e PBC + ^ + 0(l/f3 2 ) (21) 



(O)pbc = ^' FBC (22) 



Studying these corrections allows us to minimize undesired finite size effects (Ref{L2"D. 

Since the FBC partition function Eq.|| includes all possible BC, it is possible (at 
least analytically) to obtain (0) PBC directly from the FBC partition function: 

' OS(A-O) ] 

'M -<*>)„*, 

where 5(x) is the Dirac distribution. 

For most observables, (0) FBC and (0) PBC differ solely by corrections to scaling. 
This is, however, not true for the spinwave stiffness. 

Evaluating the spin rigidity using FBC is more involved: this quantity is a (spin) 
current-current corelation function and measures the phase coherence of the system. 
The standard definition of 7, namely the response of the system to a shift at the 
boundary, trivially leads to a zero value. This follows from the very implementation 
of FBC, since for any imposed shift, the dynamical variable A can adapt itself to 
absorb the shift at no cost in free energy. We need another way to compute the 
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spinwave stiffness. For instance, from Eq.[l8] the spin rigidity can be computed via 
the susceptibility of A: one should keep in mind that this expression is only valid 
for /?7 n >> 1, ...,/?7 DD >> 1, because we have extended the domain of integration 
of A' to [—00, +00]. Near a commensurate-incommensurate boundary, at least one 
of the 7 aa — > 0, and so one needs to use another procedure. 

Another limitation of this approach concerns inhomogeneous states: in that 
case one cannot use a single, constant A; for domains separated by domain walls, 
for instance, each domain wall constitutes a boundary and A depends on spatial 
coordinates. 

To overcome these problems, we will now introduce A— histograms, in conjunc- 
tion with FBC. 

B. A— Histograms 

In part |I A| we showed that the partition function with FBC is a sum over par- 
tition functions Z(A). A practical way to perform this sum is to count the number 
of configurations obtained for each of the allowed values of Ai, Ad. Since these 
quantities are defined modulo ^ , this can be easily done by histograms in A 1; A D 
which we call A— histograms. 

From the standpoint of a Monte Carlo simulation, histograms are generated as 
follows: we divide the range of variation of Ai, Ao into smaller sub- intervals. For 
each of these we store the total number of configurations n(A) having Ai,...,A£> 
within the given interval, and also the average of relevant quantities (such as the 
energy) over these configurations. 



10 



This yields the probability distribution P(A) for A, and averages O hlst (A) of 
various observables. 
-P(A) is given by 

P(K) = =$1 (23) 

where A/" is the total number of generated configurations. Similarly 
O fes *(A) is given by: 

O hist (A) = -ir£ 0(C(A)) (24) 

™( A ) C(A) 

where 2c(a) * s ^ ne sum over a ^ configurations having their boundary phase shift 
Ai, in the same given histogram interval. 

According to Eq|| P(A) in Eq.23 is proportional to e~ /3LD -^ A - ) . Thus we have 
direct access to the free energy dependence on A. A minimum of the free energy 
translates into a maximum in P(A), giving both A (and thus Qo the wavevector 
of the incommensurate structure, see below) and 7. If the peak in P(A) is sharp 
enough, thermodynamic quantities are computed for the most probable value of A; 
deviations are expected when 7 — > or when P(A) has a multi-peak structure. 



C. Thermodynamics 

From our previous discussion, we can extract relevant quantities from a Monte 
Carlo study either by using FBC without histograms, or by using FBC with his- 
tograms. This latter approach will give better results when the system is close to a 
phase transition. We can implement histogram techniques in our simulation in two 

ways: 
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1. Taking numerical derivatives of the A-histogram free energy: 
The A-histogram free energy density is obtained from: 



/(A) = -— - In (P(A)) + Constant (25) 



The zeroes of the first derivative of the free energy yield the value of A . The 
second derivatives of the free energy computed for A = A give the components 
of the spinwave stiffness 7 by Eq.|16[ The advantage of this method lies in the 
fact that it is a fast algorithm since no observables are computed. Yet, taking 
the first and second numerical derivatives of the free energy induces additional 
sources of errors in the results; however this can easily be corrected by using 
polynomial approximations in the vicinity of A . This method is well suited 
for a scaling analysis. 

2. Fluctuation-dissipation-like theorem with A-histograms: 

This approach is similar to the previous method, except that the explicit ex- 
pressions for the needed derivatives are computed analytically: using Eqs.^, |] 
we obtain Va: 



df 



d 2 f 



1 

N 



(B6 K ,) FBC (26) 



{(M*) FBO - i [(B 2 6 K ,) FBC - (B5 Kl ) 2 FBC \ } (27) 



where N is the total number of spins and: 

A = \Y1 2 J ij(tiij ■ u a ) 2 cos(y?; - + A ■ (r* - r})) (21 
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B = - Jij(uij.u a ) sin(y3i - tpj + A • (r< - r,)) (29) 

6 K , = 8(K- A') (30) 

This allows a better determination of the two derivatives compared to the 
previous method but it implies a slowing down of the algorithm, due to the 
computation of averages . 

General remarks on the A-histogram method 

• A measures the incommensurability of the system; since the partition func- 
tion is 2n/L periodic, A is obtained modulo 2n/L. Changing A to A + (^)u a 
corresponds to replacing the angular variables <pi by ipi — (^-)u a ■ r*j and leaves 
the partition function unchanged. For an homogeneous state, the helicity 
- denoted by Q - is defined modulo 2n. In the simulation we may identify its 
equilibrium value Qo with the value A corresponding to tp ~ (see Eq.^). 
Working at low T ensures that the determination of A will not drift by in- 
crements of 2it/L (a variation of A by 2n/L would entail a shift in ip by 
(2tt/L)u x , which will not occur within a reasonable number of Monte Carlo 
steps at low T). 

• According to the previous discussion, at low enough T one can construct 27r 
A— histograms: the 2tt/L periodicity does not show up in the results of the 
simulation. Closer to the transitions, energy barriers between states separated 
by 2tt/L become low and one needs to restrict the histogram to the relevant 
2n/L interval. 
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• Near phase boundaries, histogram techniques allow more accurate determina- 
tions of thermodynamic quantities. Within that framework it is also possible 
to reduce corrections to scaling, for instance by adding an extra term of the 
form //(A — A ) 2 to the partition function and thus to -P(A), and then by 
determining fi self-consistently such that corrections to scaling be minimized 
(Ref.©. 

• Using histograms can give information about the nature of the commensurate- 
incommensurate (C)-(IC) transition: in the incommensurate phase, the free 
energy displays two minima at ±A°. A first order transition will be character- 
ized by the coexistence of a third local minimum in /(A) for A = 0, at some 
characteristic temperature. 

II. APPLICATION OF THE BOUNDARY CONDITION HISTOGRAM 

TECHNIQUE 

In order to illustrate the main features of the method introduced in section |, 
we consider the row model, a frustrated anisotropic 2D XY model defined on the 
triangular latticeHi^. Only nearest neighbor sites are coupled: = —r]J (J > 0) 
for i and j along the horizontal direction and = — J otherwise. The case rj = 1 
corresponds to the fully frustrated modeli but here we will assume rj ^ 1 . At 
T = the following phase diagram is obtained. For rj < 0.5 one has a commensurate 
collinear antiferromagnetic phase. For rj > 0.5 an incommensurate spiral state is 
observed in the horizontal direction denoted by x. The entire phase diagram of the 
row model, in the (rj, T) plane, was previously studied in Monte Carlo. Since the 
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incommensurability is only present in the x direction hybrid boundary conditions 
were used: PBC in the y direction and FBC without A— histograms in the x direction. 
Three phases were identified (see Fig ([]])): a low T commensurate state (C) bounded 
by lines ALCD, a low T incommensurate state (IC) bounded by lines ALB and a 
high T paramagnetic state (P). AL is a commensurate- incommensurate transition 
line: starting from T = and increasing T at fixed rj, j xx - the x component 
of the spinwave stiffness - computed via Eq^H^, decreases, vanishes continuously 
upon reaching line AL, increases again and eventually goes back to zero on the LCD 
boundary For this particular thermodynamic path, Q x , the (incommensurate) wave 
vector in the x direction decreases and vanishes continuously on AL. 
The choice of this model stems from the fact that (%) as T approaches T^l (the 
temperature of the (C)-(IC) transition), — > and thus Eq.[T8| ceases to be valid 
(ii) for T > Tal the state of the system consists of a multidomain structure and is not 
a homogeneous (C) phase (RefflTil) and see below): (FBC), in the original formi'B, 
select A , the most probable value of A, whereas the inhomogeneous character of 
the state requires more than one A. A— histograms then prove a useful tool to study 
the thermodynamics of the row model in the vicinity of the (C)-(IC) transition. 

In the present simulation, we used A— histograms. A standard Metropolis algo- 
rithm was applied to the spin angles and to the boundary shift in the x direction. 
We worked with a 48 2 lattice and the number of MCS/spin was of order 10 5 — 10 6 . 
Typically the first 10 4 steps were discarded for equilibration. We monitored Q x , j xx , 
7 OT (the y component of the spinwave stiffness) and the staggered chiralitieJii. We 
chose f] = 0.55 and varied the temperature. A detailed analysis of the simulation is 
presented in (Ref fllif )). 
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1. At low temperature, in the incommensurate phase (IC) (see Fig ([!])), we can 
construct 27r A— histograms as discussed in the general remarks of section |. 
The corresponding histogram for Q x is shown in Fig (Q). At a given T, P(Q X ) 
displays two well defined peaks for Q x = +Q (T) and Q x = —Q (T) (the figure 
only shows positive values of Q x ). As the temperature brings the system closer 
to line AL, the peaks at Q x = +Q (T) and Q x = —Q (T) broaden and tend 
to merge. In that regime, A drifts easily and if one still uses 27r histograms 
one finds that Q x = for T > Tal (Fig (|3P). On the other hand if the 2s 
periodicity of A x is enforced, one finds that P(Q X ) shows three peaks as a 
function of Q x , for T > Tal (Figs (^) and (|5|)). The relative weight of the side 
peaks compared to that of the central peak (Q x = 0) is roughly one: 25% of 
the system is in a spiral state with wavevector Q x = +Q (T), 25% is in a spiral 
state with wavevector Q x = —Qq(T) and 50% in a collinear state. The multi- 
peak structure evokes the hysteresis often present in first order transitions but 
in fact it is due to the existence of an inhomogeneous structure for T > Tal 
(see below). 

2. At low temperature, in the incommensurate phase (IC), the determinations of 
7 XX via Eq.|18|, via 2n A— histograms or via 2s A— histograms all agree (Figs (|6]) 
and (0)). However for T > 0.10J the three curves move apart. In particular, 
for T > Tal ~ 0.19J, j xx as obtained from 2s A— histograms is zero within 
our error bars, whereas it is positive with the other methods. 

7 ra does not show any singularity : it decreases as T increases and vanishes 
on the paramagnetic boundary (line LCD of Fig (fj)). 
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For T — > Tal from above, spins undergo critical fluctuations. At a given T the 
MC process generates fluctuations in A by large increments of 2ir/L. If we use 
27r A— histograms instead of 2ir/L A— histograms (as one should) these fluctuations 
smear out the structure of the 2ir/L A— histogram (Fig([|)). The 2n A— histogram 
(Fig(|3|)) yields a single peak at Q x = 0: this value is the mean of Q x averaged over 
the 27r/L A— histogram; similarly the width of the 27T A— histogram is of the order 
of the average of Q\ over the 2ir/L A— histogram. 

The multi-peak structure of the 2tt/L A— histogram also signals the breakdown 
of the fluctuation-dissipation theorem: computing ^ for the most probable value 
of Q x (i.e Q x = 0) gives a larger result than that obtained via Eq.^7[ The former 
value of 7 XX is comparable to what one gets via 2n A— histograms. This explains 
Fig(|7p. It also explains why ^ xx determined via the most probable value of Q x is 
non zero for T > Tal- 

The fact that j xx vanishes continuously at T = Tal suggests that the 
commensurate-incommensurate transition is not first order. The structure of P(Q X ) 
and the temperature dependence of j xx for T > Tal signal in fact an inhomogeneous 
equilibrium state with striped domains in which Q x = +Qo(T), Q x = —Qq{T) or 
Q x = oi&@. We dub this phase a smectic-like phased!: in this state, the elasticity 
in the x direction vanishes, whereas it is finite in the y direction. The hamiltonian 
describing long wavelength fluctuations is smectic-likeS, i.e with a momentum de- 
pendence of the inverse propagator varying as Q\ for the x part and quadratically 
with Q y for the y part. The structure is exhibited in a snapshot of the chiralities 
of the system (Fig (||)): chiralities are positive for domains having Q x = +Qo(T), 
negative for domains having Q x = —Qo(T) and zero for walls (Q x = 0). At the 
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transition, for T = T^l, the stripe structure is stabilizedBta. 

The same characteristics are observed in the entire region ALC of the phase 
diagram (see Fig (0)). 

So, in contrast to the previous MC study, we find here four phases in the (rj, T) 
plane. In addition to the (C), (IC) and (P) phases, an inhomogeneous thermody- 
namic state is observed between the (C) and (IC) phases, bounded by lines ALC 
(see Fig ([!])). Its presence was revealed by our method. 

III. CONCLUSION 

We have presented a technique allowing to study numerically the ground state 
and the critical regime of finite size inhomogeneous or incommensurate structures. 
We illustrate its features on the commensurate-incommensurate transition of classi- 
cal 2D XY spins on the triangular lattice 
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FIGURES 





FIG. 2. Unnormalized 2tr A— histogram P{Q X ) for the wavevector Q x at various 
temperatures. Only positive values of Q x are shown, since P(Q X ) = P{—Qx)- 
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FIG. 3. Unnormalized 2tt A— histogram P{Q X ) for the wavevector Q x at various 
temperatures. For T > Tal ~ 0.19J, P{Q X ) has a single maximum at Q x = 
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FIG. 4. Unnormalized A-histo gram P(Q X ) versus Q x for rj — 0.55 and T — 0.19 J. 
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FIG. 5. First derivative of the free energy density versus Q x for r] = 0.55 and 
T = 0.19J. 
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FIG. 6. Comparison of MC data for ^y xx using A-histograms modulo 2tt or using the 
fluctuation-dissipation theorem (Eqjl8|). 
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FIG. 7. Comparison of MC data for j xx using A-histograms modulo 2n or modulo 
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FIG. 8. Snapshot of chiralities on each plaquette of a 36 2 triangular lattice, for 
T = 0.4J. Filled circles represent plaquettes with the correct sign, i.e in the same chiral 
state as at T = 0. Open circles correspond to plaquettes with the wrong sign, that is 
such that the chirality has changed compared to T = 0. Plaquettes with zero chirality 
(no symbol) are obtained in-between the two. One clearly sees a stripe structure of filled 
circles and open circles separated by domain walls of zero chirality. The horizontal (x) 
axis corresponds to the direction of the r\ bonds. 
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